% !Mode:: "TeX:UTF-8" 

\BiChapter{未来拓扑概率精确预测方法}{Exact probability computation for future connectivity}

\BiSection{引言}{Introduction}

多智能体系统的相互通信和观测是实现内部信息交流和任务协同的关键环节。%
相邻节点之间是否存在连通依赖于节点之间的相对空间位姿关系以及对应设备的通信、观测模型。%
然而多智能体系统必须保持动态响应能力以因应不同的任务需求，且智能体移动平台对应设备的相对通信和观测能力受到现实条件约束(例如平台功率限制、环境障碍遮挡等)。%
因此实际任务中，多智能体系统内部的相互通信和观测拓扑并非固定不变，而时刻变化的。


在传统的协同控制和定位任务中，系统的拓扑结构往往是当系统执行完通信和观测过程后，从通信和观测数据中提取而来。%
即当获取到的数据中存在与某一节点的通信和观测记录时，则认定该连通是存在的。%
因此在已有任务中，系统内部的观测和通信拓扑通常被视为确定且已知的。

然而系统在未来时刻的相对通信和观测拓扑在当前时刻是未知的。%
在复杂拒止环境下，协同任务的执行与完成仅分析过去和现在时刻系统和环境的状态是不够的，往往还需要对系统在未来时刻的行为进行预测，以帮助多智能体系统优化决策，从而获取最高的任务效益。%
与当前和过去时刻的拓扑结构不同的是，预测过程无法真正地执行通信和观测行为，不存在真实的通信和观测数据。因此在决策时刻，多智能体系统无法用传统方法提取未来的通信和观测拓扑。

另外，在拒止环境下，基于当前时刻的已知信息，系统未来时刻通信和观测拓扑的分布本质上是随机的。%
由于拒止环境下智能体无法获取准确的定位信息，因此必须通过融合对环境和自身状态的测量数据来估计其自身位姿信息。%
一方面传感器的测量过程并不精确，而是受到测量噪声的干扰；另一方面，多智能体内部也存在未建模的动力学扰动。%
节点对自身真实位置的估计通常是一个概率分布函数，需要一个位姿状态均值和分布协方差矩阵来表示。%
虽然过去和当前时刻的通信和测量拓扑可以通过实际已经获得的通信和观测记录来确定，但是系统对未来时刻通信和观测拓扑的预测却无法如此获得。%
任意两个节点之间的连通存在与否，是与两个节点未来位姿的预测分布以及通信和观测设备的模型相关的一个随机事件。

目前已有文献对预测未来时刻系统的通信和观测拓扑的方法主要是采用最大似然(Maximum Likelyhood)
假设，即仅考虑具有最大可能性或最大概率的未来表现。%
当系统的干扰模型全部为高斯白噪声时，任意两个节点之间的连通将仅依赖于未来时刻节点预测分布的均值以及通信和观测设备的模型。%
不过该方法仅考虑了一种确定的未来拓扑，无法解释系统所有可能的未来行为。

本章研究在当前时刻根据已知信息计算任意节点之间在未来时刻的连通概率。%
与最大似然假设不同，本章不仅仅使用节点在未来时刻的预测位姿均值，而是考虑两个节点完整的预测位姿分布（包括均值和分布的协方差矩阵）。%
假设系统的干扰均为高斯模型，则针对一类常见的通信测量模型——圆盘模型，任意两个节点在未来时刻的连通存在与否，不再是最大似然假设中或存在或不存在的伯努利分布，而是依特定概率存在的正太二次型分布。本章给出具体的建模过程，并基于正太二次型多项式展开定理，提出一类精确计算未来连通度概率的方法，为后文基于概率拓扑的主动定位规划问题提供研究基础。

\begin{figure}[ht]
	\centering
	\includegraphics[width = \textwidth]{figures_PROB/Task}
	\bicaption[fig. PROB_task]{}{概率连通问题示意图。}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

\BiSection{问题建模}{Problem formulation}

本章假设多智能体系统由$N$个智能体和$O$个固定锚节点组成。其中领导者为可移动的智能体，均搭载了相同的相对测量和通信设备。它们通过执行相互通信和相对观测行为进行组网，共同组成一个多智能体系统，在GNSS拒止的环境中执行相关任务。因此，智能体无法获取自身的精确位置，需要根据配备的传感器数据(如里程计、相对距离测量等)以及通信设备中解码的其他智能体传输的数据对自身位置进行估计。锚节点为一类已知自身精确位置信息的节点，它们装备有与智能体相同的测量和通信设备，可以对其有效覆盖区域内的其他节点进行相对测量并持续广播自身位置信息。用$\vect{p}_i^{k}$ 和 $\vect{W}_j$代表第$i$个智能体在第$k$时刻的真实位置和第$j$个锚节点的真实位置。

{\textbf{图论：}}使用一个时变无向拓扑图$\mathcal{G} = \{\mathcal{V}, \mathcal{E}(k)\}$描述多智能体系统中的相对观测和通信关系。图中的节点包括智能体和锚节点，统一使用$v_i$或数字索引$i$表示，因此$\mathcal{V}:=\{ v_i| \forall i \in \mathbb{R}^+, i \le N+O \}$为所有节点的集合，包括智能体集合$\mathcal{V}_R$，以及锚节点集合$\mathcal{V}_A$。$\mathcal{E}(k)$为两个相邻节点之间存在的相对测量和通信连通的集合。根据式~\ref{eq. PROB_regualtionFunction01}的定义，$\mathcal{E}(k)=\{ v_i \times v_j | \forall  s_{i,j}(k)=1, i \neq j \}$,式中$s_{i,j}(k)$为$s_{\rho}(d_{i,j})$的简化表述。当且仅当$s_{i,j}(k)=1$时，两个节点之间互为邻接节点。使用$\vect{N}_i(k)=\{ v_j | \forall s_{i,j}(k) \in \mathcal{E}(k), s_{i,j}(k) = 1\}$为节点$v_i$在$k$时刻可以进行相互测量和通信的邻接节点集合,用$n_i^k$表示$k$时刻第$i$个智能体邻接节点的个数。

智能体的运动模型为：
\begin{equation} \label{eq. PROB_motionModel}
	\vect{p}_i(k+1) = f_i \left[ \vect{p}_i(k),\vect{u}_i(k),\vect{w}_i \right],
\end{equation}
式中$\vect{p}_i(k)$为第$i$个智能体在$k$时刻的位置状态，$\vect{u}_i(k)$表示第$i$个智能体在$k$时刻的控制输入，$\vect{w}_i$为该智能体运动过程受到的外部随机扰动。假设该随机扰动服从已知的正太分布$\vect{w}_i \sim \mathcal{N}(\vect{0}, \vect{R}_i)$。

假如在$k$时刻，第$i$个智能体和第$j$个智能体或锚节点之间存在相互测量，则观测数据$\vect{\vec{z}}_{ij}^k$为：
\begin{equation} \label{eq. PROB_observationModel}
	\vect{\vec{z}}_{i,j}(k) = h_{i,j} \left[ \vect{p}_i(k),\vect{p}_j(k), \vect{v}_{i,j} \right],
\end{equation}
式中$\vect{v}_{i,j}$为相对测量设备中存在的随机观测误差，假设该误差依然服从一个已知的正太分布$\vect{v}_{i,j} \sim \mathcal{N}(\vect{0}, \vect{Q}_i)$。
观测量$\vect{\vec{z}}_{ij}$的具体数据格式与智能体携带的相对观测设备的种类有关。例如当相对观测设备为声纳、激光测距等设备时，仅测量两个节点之间的相对距离，则$\vect{\vec{z}}_{ij}$仅为标量。当相对观测设备为雷达等即可测量相对距离又可测量相对角度的设备时，$\vect{\vec{z}}_{i,j} \in \mathbb{R}^{2}$为二维向量。后文将混合使用上标$(\cdot)^k$和括号$(\cdot)(k)$两种符号来表示第$k$个时刻。

仅依赖于式~\eqref{eq. PROB_observationModel}无法精确描述任意时间两个节点之间由相对观测过程而建立的联系。下面将节点$i$在$k$时刻对任意节点$j$的观测行为定义为事件。由于存在如下两点原因，

1）：节点的运动模型和观测设备的观测模型都存在误差；

2）：观测设备的观测能力受到实际因素地限制，节点之间的观测连通无法保持全连接、也无法维持固定的连通关系，而是会根据节点之间的相对位置关系而时刻改变；

因此观测事件本质上是一组包含三个变量的元胞，即$\vect{z}_{ij}(k) = < \vect{\vec{z}}_{ij}^k, s_{ij}^k, p(s_{ij}^k) >$，其中$s_{ij}^k$为图论中定义的连通度二进制变量，当$s_{ij}^k=1$时，表示$i$和$j$两个节点之间存在观测连通，否则$s_{ij}^k = 0$；$p(s_{ij}^k)$表示连通度存在与否的概率；$\vect{\vec{z}}_{ij}^k$仅在$s_{ij}^k=1$时有效，代表当两个节点$i,j$存在观测连通时，观测设备实际获得的观测数据。

假设集群已经执行了$k$次控制输入$\vect{U}^{0:k-1}$到达了$k$时刻，并且集群在每一时刻都会进行相对测量操作，得到观测数据$\vect{Z}^{0:k}$。则对于过去时刻的观测$\forall t \le k-1$，每一组观测元胞中连通变量$s_{ij}^t$和连通概率$p(s_{ij}^t)$的值都可以从观测数据$\vect{Z}^{0:k}$中提取得到，即
\begin{equation}
	\left\lbrace
	\begin{array}{ll}
		s_{ij}^t=1,p(s_{ij}^t=1)=1,p(s_{ij}^t=0)=0 , \quad &  \text{如果$\vect{z}_{ij}^t \in \vect{Z}^{0:k-1}$}\\
		s_{ij}^t=0,p(s_{ij}^t=0)=1,p(s_{ij}^t=1)=0, \quad & \text{否则}, 		
	\end{array}			
	\right.
\end{equation}
然而，如果给定一组对未来$L$时刻的控制输入$\vect{U}^{k:k+L}$，则在该控制输入作用下，未来的运动过程和相对观测过程都尚未发生，因此未来时刻每一个观测事件$\vect{{z}}_{ij}^t, \forall k+1 \le t \le L$ 都是未知的，因此节点无法再依照上式的判定准则预测未来时刻的连通、连通概率以及实际观测值。

本章重点研究如何基于$k$时刻的已知数据，例如节点的估计位置和估计不确定性，以及规划的控制输入等，预测任意节点之间观测事件的连通变量，以及连通概率。本文考虑一类典型的圆盘模型模拟节点之间的观测连通变量，
即假设智能体携带的通信和相对测量设备因受到实际约束而存在最大工作半径，并且假设智能体和锚节点上的测量设备和通信设备的工作半径是相同的，全都以$\rho$表示。则该模型的数学表达式为：
\begin{equation}	\label{eq. PROB_regualtionFunction01}
	s_{\rho}(d_{ij}(k)) = \left\lbrace
	\begin{array}{lr}
		1, & d_{ij}(k) \le \rho \\
		0, & d_{ij}(k) > \rho, 		
	\end{array}			
	\right.
\end{equation}
式中$d_{i,j}(k)$为两个节点之间真实的相对距离：
\begin{equation*}
	d_{i,j}(k) = \left\lVert \vect{p}_i(k) - \vect{p}_j(k) \right\rVert_2.
\end{equation*}

下面给出后文将使用的一个重要定理\cite{provost1992quadratic}，

\begin{theorem}  \label{th.PROB_PowerSeriesExpansion}
	给定一个$p$维的多元正太分布向量$\vect{X} \sim \mathcal{N}(\vect{\mu}, \vect{\Sigma})$, $\vect{\Sigma} > 0$，定义它的二次型$Y = \vect{X}^T\vect{A}\vect{X}$，其中矩阵$\vect{A}$式一个正定对称矩阵，即$\vect{A} = \vect{A}^T > 0$。用$\vect{\lambda} = [\lambda_1,...,\lambda_p]^T$表示矩阵$\vect{\Sigma}^{\frac{1}{2}}A\vect{\Sigma}^{\frac{1}{2}}$的特征值，且特征向量矩阵为$\vect{P},\vect{P}^T\vect{P} = \vect{I}$， 即 $P^T\vect{\Sigma}^{\frac{1}{2}}A\vect{\Sigma}^{\frac{1}{2}}P = diag(\vect{\lambda})$。则定义向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$以及给定一个正数$y$，则二次型随机变量$Y$的概率函数可以通过如下的多项式序列展开进行计算：
	\begin{equation} \label{eq.PROB CDP_infinity}
		p(Y\le y) = \sum_{w=0}^{\infty}(-1)^w c_w \frac{y^{\frac{p}{2} + w}}{\Gamma(\frac{p}{2}+w+1)}, \quad 0<y<\infty,
	\end{equation}
	式中的参数序列$c_w$的计算方法为：
	\begin{equation} \label{eq.PROB_c_w}
		\left\lbrace
		\begin{aligned}
			c_0 &= exp(-\frac{1}{2}\sum_{j=1}^{p}b_j^2)\prod_{j=1}^p(2\lambda_j)^{-\frac{1}{2}}, \\
			c_w &= \frac{1}{w}\sum_{r=0}^{w-1}d_{w-r}c_r, \quad w \ge 1,
		\end{aligned}
		\right.
	\end{equation}
	以及
	\begin{equation}\label{eq.PROB_d_w}
		d_w = \frac{1}{2} \sum_{j=1}^p (1-wb_j^2)(2\lambda_j)^{-w},\quad w \ge 1.
	\end{equation}
	式中，$\Gamma(\cdot)$函数为伽马函数，当输入为正整数时，如$n\in \mathbb{R}^+$，其输出等于$(n-1)!$；$b_j$表示向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$中的第$j$个元素。
\end{theorem}

\BiSection{基于有限项近似的概率计算方法}{arg2}
本节设计一种多智能体系统在获取预测状态后两个节点在未来时刻存在通信和观测连通的概率计算方法。该方法首先基于通信和观测设备的圆盘模型，将连个节点之间的连通转变为一个正定二次型随机变量的概率分布，然后使用定理~\ref{th.PROB_PowerSeriesExpansion}计算未来时刻存在连通的概率。由于定理~\ref{th.PROB_PowerSeriesExpansion}中概率的计算公式需要累加无穷多个多项式，因此为了易于计算，仅选用有限个多项式的累加来近似真实概率。然后为了减少因有限项近似引入的截断误差，本节提出了三个计算过程的修正方法。最后进行数值仿真，验证本方法的有效性。

\BiSubsection{有限项近似}{Finite-term Approximation} 
为了在$k$时刻评估控制输入$\vect{U}^k$，多智能体系统在调用合适的估计引擎后，得出了对下一时刻，即$k+1$时刻，系统状态的预测均值向量$\vect{\tilde{X}}^{k+1}$以及协方差矩阵$\vect{\Sigma_{\tilde{X}^{k+1}}^k}$。此时为了预测在$k+1$时刻任意两个节点，例如节点$\vect{p}_1^{k+1}$和节点$\vect{p}_2^{k+1}$，之间通信和观测连通存在的概率，可以从预测值中提取对应节点在$k+1$时刻的预测状态分布。由于本章默认采用高斯噪声假设，因此这两个节点在未来时刻的预测状态也是两个高斯分布函数，即$\vect{p}_1^{k+1}\sim \mathcal{N}(\vect{\tilde{p}}_1^{k+1}, \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}})$以及$\vect{p}_2^{k+1}\sim \mathcal{N}(\vect{\tilde{p}}_2^{k+1}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}})$。因此通过减法定义的上述两个正太分布的差也将符合高斯正太分布函数，即$\Delta \vect{p}_{12}^{k+1} = [\left( \Delta \vect{p}_{12}^{k+1} \right)_x, \left( \Delta \vect{p}_{12}^{k+1} \right)_y]^T \sim \mathcal{N}(\vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}})$。

定义节点$\vect{p}_1^{k+1}$和节点$\vect{p}_2^{k+1}$在未来的$k+1$时刻的真实距离为$$d_{1,2}(k+1) = ||\vect{p}_1^{k+1} - \vect{p}_2^{k+1}||_2,$$
则将该距离视作随机变量，它的平方可以表述差的正太分布$\Delta \vect{p}_{12}^{k+1}$的二次型的形式：
$$Y = d_{12}^2(k+1) = \left( \Delta \vect{p}_{12}^{k+1} \right)_x^2 + \left( \Delta \vect{p}_{12}^{k+1} \right)_y^2 = \left( \Delta \vect{p}_{12}^{k+1} \right)^T \vect{A} \Delta \vect{p}_{12}^{k+1},$$
式中矩阵$\vect{A} = \vect{I}$为单位阵。

基于通信设备和相对测量设备的圆盘模型，两个节点在未来的$k+1$时刻存在通信和观测的概率，即$s_{12}^{k+1} = 1$的概率，可以等价地转换为
\begin{equation*}
	p\left( s_{1,2}^{k+1} = 1 \right) = p\left( Y\le \rho ^2\right).
\end{equation*}
式中，$\rho$为通信设备和相互测量设备最大的工作半径。

因此对未来时刻连通概率的预测可以转换为两个节点相对距离平方的概率分布，即可以使用定理~\ref{th.PROB_PowerSeriesExpansion}来计算相应的概率。当考虑二维空间的任务规划时，所有节点的状态均为二维向量，因此定理~\ref{th.PROB_PowerSeriesExpansion}中的维度$p=2$。又由于该概率计算公式需要累加无穷多个多项式，在数值计算中是不可行的，因此本章提出仅使用有限项，即前$w = {1,2,...,w_F}$项来做近似计算。由此式~\eqref{eq.PROB CDP_infinity}可变为：
\begin{equation} \label{eq. PROB_CDP_finite}
	p(s_{12}^{k+1} = 1) = p\left( Y \le \rho ^2\right) =  \sum_{w=0}^{w_F}(-1)^w c_w \frac{ （(\rho^2)^{w + 1}}{(w+1)!},
\end{equation}
然而这种通过截取有限项来近似地计算概率的方法必然会引入计算误差。因此上式在应用上的难点是如何保持算法在计算上的可行性的同时减少由有限项近似引入的截断误差。

\BiSubsection{截断误差分析与修正方法}{arg2}

本节针对有限项近似对概率计算公式\eqref{eq.PROB CDP_infinity}的影响进行讨论，通过仿真实验结果观察到两个主要的问题所在，并根据这些问题提出了三种修正方法。

\begin{figure}[ht]
	\centering
	\includegraphics[width = \textwidth]{figures_PROB/PerformanceEvaluationTwoNodes1}
	\bicaption[fig. PROB_PerEvaTwoNodes]{}{概率连通问题示意图。}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

如图~\ref{fig. PROB_PerEvaTwoNodes}所示，本节使用有限项近似的计算结果与基于正太分布抽样点的频率对比，用以侦测有限项近似对频率计算公式~\ref{eq.PROB CDP_infinity}的影响。此处选取两个随机分布节点相对观测的场景。场景一的两个节点正太分布如图~\ref{fig. PROB_PerEvaTwoNodes} A 所示，场景二如图\ref{fig. PROB_PerEvaTwoNodes} D 所示。两个场景中，节点之间的平均距离均为$5m$。两个节点之间的相对距离的分布如图 B 和 E 所示。两个场景的区别在于场景二两个节点之间相对距离的分布比场景一更为稠密。图 C 和 F 为两个节点之间的相对距离的密度分布示意图，两种仿真场景总共抽样了10000个样本，纵坐标为两个节点相对距离与横坐标相等的样本数。 图 G 和图 H 为两个场景中，当观测和通信半径$\rho$从0变化到10的时候，两个节点之间的通信和观测连通存在的真实概率（以10000次抽样的频率为真值参考），以及使用不同$w_F$的有限项近似公式得到的近似概率。

从图~\ref{fig. PROB_PerEvaTwoNodes} G-H中可以观察到两个主要的误差来源。首先是观察子图 E 中当$w_F = 50$和 $w_F = 100$时，有限项近似公式计算得到的近似概率会随着通信和观测设备工作半径$\rho$的增加而产生跳变（$w_F = 50$的变化曲线在约$\rho>6$时突然增加超出概率的$[0,1]$域值，$w_F = 100$的变化曲线在$\rho>8$时突变为负值）。有限项近似的计算方法仅可稳定且准确地计算工作半径$\rho$小于突变点的概率值。突变的方向与使用的多项式个数$w_F$相关。通过观察式~\eqref{eq.PROB CDP_infinity}，该误差的产生源自累加项中存在的$(-1)^k$。该跳变数列的存在使得两个相邻的累加项随设备半径$\rho$的变大而趋向于完全相反增长方向。因此当累加无穷多个多项式计算时，得出的概率可以被限制在合理的$[0 1]$之间。然而当选取有限个多项式的累加作近似计算时，由于忽略了后续的多项式，因此导致朝向一个方向增长的计算数值得不到平衡，因此整体的概率计算向该方向无限增长，从而超出合理范围。

有限项近似导致的第二个计算误差如图~\ref{fig. PROB_PerEvaTwoNodes} H 所示，即当相对距离的分布较为集中，从而导致距离分布的协方差矩阵较小或者近似奇异时，有限项近似方法在临界半径范围内（图 H 中约 $4 < \rho < 6$），无论如何选取累加的多项式数目$w_F$，均无法估计连通存在的概率。
该问题的由来可追溯至定理~\ref{th.PROB_PowerSeriesExpansion}中参数序列$c_w$和向量$\vect{b}$的定义。根据前一节的建模过程以及定理~\ref{th.PROB_PowerSeriesExpansion}，$c_0$的计算需要用$\vect{\Sigma}^{\frac{1}{2}}\vect{A}\vect{\Sigma}^{\frac{1}{2}}$的特征根的倒数。而当二次型随机变量为节点之间相对距离的平方时，$\vect{A} = \vect{I}$，因此$c_0$需要用到$\vect{\Sigma}$特征根的倒数。
另外，向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$也需要用到协方差矩阵$\vect{\Sigma}$的逆。因此当该协方差矩阵的特征值较小或者矩阵近似奇异时，参数序列$c_w$的向量$\vect{b}$的值将被异常放大，从而导致有限项近似的计算偏离正确概率。

\begin{figure}[ht]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/3_delta_figure}
	\bicaption[fig. PROB_3delta]{}{$3\delta$修正方法示意图}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

针对以上两种计算过程引入的误差，本节提出三种修正方法。
\begin{itemize}
	\item[1)]  {\hei $3 \delta$区间校正}: 基于高斯噪声假设，两个节点的分布$\vect{p}_1^{k+1}$和$\vect{p}_2^{k+1}$均为正太分布，因此两个节点之间通过减法定义的矢量$\Delta \vect{p}_{12}^{k+1}$也是一个二维正太分布的向量。如图根据高斯正太分布函数的性质可知，$\Delta \vect{p}_{12}^{k+1}$分布中$99.7\%$的抽样点均将落入其$3\delta$置信区间中。因此，此处的$3 \delta$区间校正方法是指在获取节点之间的减法矢量$\Delta \vect{p}_{12}^{k+1}$的分布后，通过提前计算该分布的$3\delta$区间，得到边界$[\rho_1,\rho_2]$，则当节点的通信观测设备的工作半径$\rho<\rho_1$时，将节点之间存在连通的概率直接赋值为0，否则当$\rho>\rho_2$时，将概率直接赋值为1。

	\item[2)] {\hei 突变校正}：经过$3 \delta$区间校正后，第一类突变误差中向正无穷和负无穷方向的突变都将被直接赋值为0或1。然而当$\rho$继续增加时，仍可能会出现从1直接跳变为0的情况（如图~\ref{fig. PROB_PerEvaTwoNodes} H 中$w_F=100$的变化曲线）。观察到此类突变总是发生在工作半径$\rho$逐渐增加到比较大的数值后，因此突变校正的方法主要是当通信和观测设备的工作$\rho$大于两个节点之间的平均半径后，通过有限项近似计算出的节点之间存在连通的概率却为零或负值，此时则认为产生了此类突变，因此直接将连通概率设定为1。
	
	\item[3)] {\hei 协方差矩阵近似膨胀}：以上的两种校正方法均是针对突变误差采取的修正措施。在多智能体系统的协同任务中，处于临界状态的连通往往决定了系统的最佳性能表现。因此关注临界的通信和工作半径状态下节点之间连通的计算过程是十分重要的。由以上分析可知，有限项近似的概率计算方法在相对距离分布的协方差矩阵较小或者近似奇异时表现不佳。因此，解决该问题的思路是将较小或者近似奇异的矩阵进行等比例膨胀，将膨胀后正常的矩阵带入有限项近似的概率计算公式，作近似计算。
\end{itemize}

\begin{figure}[ht]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/ACE}
	\bicaption[fig. PROB_ACE]{}{协方差近似膨胀校正方法}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

下面将详细介绍协方差近似膨胀方法的设计过程。如图~\ref{fig. PROB_ACE}所示，红色实线的椭圆$E_{r1}$代表两个节点相对距离分布协方差矩阵$\vect{\Sigma}$的$3\delta$置信域。此时两个节点上配备的通信和测量设备的工作半径为$\rho_1$，则以坐标原点为中心，$\rho_1$为半径的通信测量圆与$E_{r1}$相交于$D$和$F$点。假设$E_{r1}$内部存在一个小的椭圆$E_{b1}$与半径为$\rho_1$通信测量圆相切于$E$点。由于$E_{r1}$所代表的协方差矩阵较小或者是近似奇异，因此将$E_{r1}$在$O_1$点等比例放大$\beta$倍，获得新的置信域$E_{r2}$。通过等比例放大，内部的相切椭圆$E_{b1}$也被放大为$E_{b2}$。$E_{r2}$相比于$E_{r1}$，其所对应的协方差矩阵经过成倍扩大后，不会存在奇异现象。假设存在一个新通信观测圆，以$\rho_2$为工作半径，与$E_{r2}$相交于$A,C$两点，与$E_{b2}$相切于$B$点。

则由相关定义可知，两个节点之间存在通信和观测连通的概率等于节点之间相对距离小于设备工作半径，即图中通信观测圆与置信椭圆之间的交集部分，其数学表述为：
 \begin{equation*}
	p(s_{12}^{k+1} = 1) \approx \frac{S(\widearc{DEF} \cap E_{r1})}{S(E_{r1})} 
\end{equation*} 
式中$S(\cdot)$为计算输入参数的面积。为了将合适的膨胀后的协方差带入有限项近似的概率计算公式中，需要寻找一个合适的$\rho_2$使得：
 \begin{equation*}
	\frac{S(\widearc{DEF} \cap E_{r1})}{S(E_{r1})}  \approx \frac{S(\widearc{ABC} \cap E_{r2})}{S(E_{r2})}
\end{equation*} 
本节提出一种简易方法来近似地计算一个合适的$\rho_2$。该方法假设三组点群：$\{ O_1, F, C \}$,$\{ O_1, D, A \}$和$\{ O_1, E, B, O \}$中每一组里的点都是相互共线的。则根据几何相似关系，可得
\begin{equation} \label{eq. PROB_NewRadius}
	\rho_{2} \approx (1 - \sqrt{\beta})d_m + \sqrt{\beta}\rho_{1}.
\end{equation}
式中$d_m$为两个节点之间的平均距离。此过程的修正方法中膨胀系数$\beta$的值实际上可以随意指定。后文确定$\beta$的方法是确保膨胀后的矩阵最小特征根大于1。

\BiSubsection{修正后的有限项近似概率算法}{arg2}
基于前一节的分析，本节给出融合上述三项修正后的有限项近似的概率计算算法~\ref{alg. finite_term_approximation}。该算法的输入为两个节点对其在$k+1$时刻位置的估计分布(均值$\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1}$以及协方差矩阵$\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$)，以及节点的通信和观测设备最大工作半径$\rho_{1}$，选定的多项式累加项的数目$w_F$。输出为在$k$时刻两个节点对它们在$k+1$时刻存在通信和观测连通的概率，即$p( s_{12}^{k+1} = 1 ) $。

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwResult{ $p( s_{12}^{k+1} = 1 ) $  }
	\KwData{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}, \rho_{1}, w_F$ \;}
	
	$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
	$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
	
	$[\rho_{3\delta}^{min} \rho_{3\delta}^{max}] \leftarrow 3\delta$ 区间计算 ( $\vect{\Sigma}_{\Delta}$ )  \;
	\If{$\rho_{1} < \rho_{3\delta}^{min}$}{
		\textbf{Return} $p(s_{12}^{k+1} = 1) = 0$ \;
	} 
	\If{$\rho_{1} > \rho_{3\delta}^{max}$}{
		\textbf{Return} $p(s_{12}^{k+1} = 1) = 1$ \;
	} 
	$\beta = \dfrac{1}{min(eig(\vect{\Sigma}_{\Delta}))}$ \;
	
	\eIf{$\beta < 1$}{
		$\beta$ = 1, $\rho_{2} = \rho_{1}$ \;}
	{
		$\rho_{2} \leftarrow \text{式}~\eqref{eq. PROB_NewRadius}$ \;
	}
	$\vect{\Sigma} = \beta \times \vect{\Sigma}_{\Delta}$\;
	$\vect{b} = \vect{P}^T\vect{\Sigma}^{-\frac{1}{2}} \vect{\mu}_{\Delta};\vect{\lambda} = \vect{P}^T \vect{\Sigma} \vect{P}$, $\vect{P}^T\vect{P} = \vect{I}$\;
	
	\While{$r \le w_F$}{
		$c_r \leftarrow \text{\eqref{eq.PROB_c_w}}, r = 0,...,w_F$\;
		$d_r \leftarrow \eqref{eq.PROB_d_w}, r = 1,...,w_F$\;
		$y_r = (-1)^r \dfrac{\left( \rho_{2}^2 \right)^{r+1}}{(r+1)!}, r = 0,...,w_F $\;
	}
	$ p( s_{12}^{k+1} = 1 )  \leftarrow \eqref{eq. PROB_CDP_finite}$ \;
	\If{$p( s_{12}^{k+1} = 1 )  \le 0 $ $\&$ $\rho_{1} > ||\vect{\mu}_{\Delta}||_2$}{
		$p( s_{12}^{k+1} = 1 )  = 1$ \;
	} 
	\textbf{Return} $Pr(Y<\rho_{old}^2)$	
	\AlgoBiCaption{经过修正后的有限项近似概率计算算法}{Controller for $i$-th robot at timestep $k$}
	\label{alg. finite_term_approximation}
\end{algorithm}

该算法按照从前到后的顺序：$1-2$进行数据预处理，从预测的位置分布中得到相对距离的高斯分布函数。$3-9$行，首先计算协方差矩阵对应的$3\delta$区间，然后按照$3\delta$区间校正的过程对未来连通的概率值作出前期校正。此后的$10-16$行属于协方差近似膨胀的校正过程，主要任务是确定膨胀系数$\beta$，以及计算膨胀后的等效通信观测半径$\rho_{2}$。$17-23$行根据定理~\ref{th.PROB_PowerSeriesExpansion}以及有限项近似的计算公式~\eqref{eq. PROB_CDP_finite}进行累加计算，得到有限项近似计算的概率值。$24-26$行进行最后的突变校正过程，检测是否有突然跳变为负值概率的情况。

\BiSubsection{实验仿真验证}{Simulation results}
本节对算法~\ref{alg. finite_term_approximation}进行数值仿真验证。仿真分为对两个节点单个未来连通的预测以及同时对多个节点之间连通的预测。

首先，使用修改后的概率预测算法对两个节点单个连通的预测结果如图~\ref{fig. PROB_PerEvaTwoNodes} I-J 所示。将I-J与未修改前的有限项近似计算公式~\eqref{eq. PROB_CDP_finite}的计算结果 G-H 进行比较可以看出，算法~\ref{alg. finite_term_approximation}有效地改善了有限项近似的概率预测结果。虽然当多项式累加数目$w_F$较少时，依旧存在概率计算值跳变的情况，选取较大的$w_F$可以有效地避免此类问题。

下面测试使用算法~\ref{alg. finite_term_approximation}同时预测多个节点之间的连通。仿真场景设置为3个节点，用编号$\{ 1,2,3 \}$表示，以及一个锚节点$\{ A \}$。则任意仿真时刻，多智能体系统中最多存在六个通信和相对测量连通。设$k$时刻六条通信测量链之间的相对距离分布如表~\ref{tab. PROB_initialDistance}所示，节点上携带的通信和测量设备的最大工作半径设为5$m$。假设在$k$时刻，所有节点都可以获取其位置信息，即所有节点在$k$时刻对自身位姿估计的协方差矩阵全部为零矩阵。已知节点的运动过程受到的噪音协方差$\vect{R}_i，i=1,2,3$均为$diag([0.04,0.04])$。此处选取较小的运动过程噪声的目的是使得节点对其$k+1$时刻预测状态的协方差矩阵较小或者近似奇异，因此可以触发算法~\ref{alg. finite_term_approximation}中的协方差近似膨胀修正。
\begin{table}[htbp]
	\bicaption[tab. PROB_initialDistance]{}{多节点$k$仿真时刻的距离分布}{Table$\!$}{Distance for the simulation case of multiple connections at $k$ step}
	\vspace{0.5em}
	\centering
	\wuhao		
	\begin{tabular}{ccccccc}
		\toprule[1.5pt]
		连通编号 & 1 & 2 & 3 & 4 & 5 & 6 \\
		节点对 & \{$A$,$r_1$\} & \{$A$,$r_2$\} & \{$A$,$r_3$\} & \{$r_1$,$r_2$\} & \{$r_1$,$r_3$\} & \{$r_2$,$r_3$\} \\
		初始距离(m) & 3 & 5.83 & 4.47 & 3.61 & 4.12 & 3.16 \\
		\bottomrule[1.5pt]
	\end{tabular}
\end{table}

具体的仿真过程如下：
\begin{itemize}
	\item[1):] 首先所有的运动节点均在当前时刻（即$k$时刻）规划一组控制输入$\vect{U}^k = [ \vect{u}_1^k, \vect{u}_2^k,\vect{u}_3^k ]$,其中$abs(\vect{u}_i^k)<1, i =1,2,3$。
	\item[2):] 使用拓展卡尔曼滤波，执行先验更新过程。因为对未来行为的预测无实际的测量数据，因此节点对自身位置状态的预测仅使用$k+1$时刻的先验估计数据。然后使用修正后的有限项近似算法~\ref{alg. finite_term_approximation}分别计算六条通信观测链在下一时刻存在的概率。
	\item[3):] 从噪音协方差矩阵 $\vect{R}_i,i=1,2,3$ 中随机抽样出噪音，施加于节点的运动过程，观察六条通信观测链的存在与否。该步骤重复1000次，通过统计每个通信观测连存在的频次来模拟该连通在未来存在的概率。这个频率数据被视为该连通存在概率的真值。
	\item[4):] 将步骤2)中计算得到的频率于步骤3)中统计得到的频率做差，得出算法~\ref{alg. finite_term_approximation}对未来连通概率的预测误差。
	\item[5):] 将步骤1)-4)重复200次。
\end{itemize}
200次的预测误差分布如图~\ref{fig. PROB_ErrorDistributionMultiConnections}所示：
\begin{figure}[ht]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/PerformanceMultiConncetions}
	\bicaption[fig. PROB_ErrorDistributionMultiConnections]{}{200次仿真中多节点环境概率预测误差分布}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

观察仿真结果可以看出，首先六条通信观测链的概率预测误差上限不超过$7\%$，因此可以验证本节所提出的算法的有效性。其次，节点对之间的相对距离越接近设备工作半径的临界，算法~\ref{alg. finite_term_approximation}的预测误差越是显著。

\BiSection{基于自适应参数选取的改进概率预测计算方法}{arg2}

前一节基于有限项近似的方法，通过提出三项基于经验的修正方法，得到了一种较为稳定的连通概率预测算法。然而前一种方法存在如下三方面的问题：
\begin{itemize}
	\item[1):] 三项经验启发式的修正方法的普适性和稳定性没有保证。
	\item[2):] 有限项近似中累加的多项式数目$w_F$在算法~\ref{alg. finite_term_approximation}中是提前选定的，对协方差矩阵变化的适应性不高。
	\item[3):] 协方差近似膨胀校正过程中存在较多的近似过程，降低了算法的预测性能。
\end{itemize}

因此本节从理论上对有限项近似算法的误差根源，即截断误差，作进一步分析验证。通过严格的理论推导得出了任意的期望截断误差限制下，参数$w_F$的自适应选取方法，从而使得算法可以针对不同的相对距离分布选取不同的累加项数目。此外，本节还将对协方差近似膨胀的校正过程进行改进，提出一种基于共线平移的膨胀策略，从而避免比必要的近似过程，从而提升预测性能。

在进行主要截断误差分析之前，需要声明下面的引理，
\begin{lemma} \label{lem. PROB_factorial_lemma}
	给定一个正整数 $n$, 它的阶乘有如下关系式：
	$$ \left( \frac{n}{\mathit{e}} \right)^n < n! < \mathit{e} \left( \frac{n}{2} \right)^n$$
	式中$\mathit{e}$为自然常数。
\end{lemma} 

\BiSubsection{截断误差稳定性分析}{arg2}

通过对比定理~\ref{th.PROB_PowerSeriesExpansion}中原始的无穷项累加公式~\ref{eq.PROB CDP_infinity}和有限项近似公式~\ref{eq. PROB_CDP_finite}可以得到截断误差的表达式为：
\begin{equation} \label{eq. PROB_truncationError}
	\begin{aligned}
		\Delta p(s_{1,2}^{k+1}) = \sum_{w=w_F + 1}^{\infty}(-1)^{w} c_w \frac{(y)^{w+1}}{(w+1)!} .
	\end{aligned}
\end{equation}
式中，$y=\rho^2$。

虽然该截断误差也是无穷多个多项式的累加和，下面将通过理论推导证明，通过选取合适的$w_F$时，可以任意地减少该截断误差。首先，推论~\ref{cor. PROB_dw} 总结了公式~\eqref{eq.PROB_d_w}定义的参数序列$d_w$的性质如下：
\begin{corollary} \label{cor. PROB_dw}
	给定一个二维正太向量$\vect{X}\sim\mathcal{N}(\vect{\mu},\vect{\Sigma})$，用$\vect{\lambda} = [\lambda_1,\lambda_2]^T$表示协方差矩阵$\vect{\Sigma}$的特征值，矩阵$\vect{P}$为协方差矩阵$\vect{\Sigma}$的特征向量矩阵。定义向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{\frac{1}{2}}\vect{\mu} = [b_1,b_2]^2$，则当下列条件满足时，
	\begin{equation*}
		\left \lbrace
		\begin{aligned}
			&2\lambda_j > 1  \\ 
			&w \ge \frac{1}{b_j^2} + \frac{1}{2\lambda_j - 1} \\
		\end{aligned}	
		\right.
		,\forall j \in \{ 1,2 \}
	\end{equation*}
	公式~\eqref{eq.PROB_d_w}定义的参数序列$d_w$的极限为零，即
	$$\mathop{\lim}_{k\rightarrow \infty} d_k = 0,$$
	且序列$d_w$是有界的，即存在正的实数$\bar{d}$,使得
	$$|d_k|<\bar{d},\forall k \in \mathbb{R}^+.$$
\end{corollary}

\begin{proof}
	该推论的证明参见附录~\ref{appendix. PROB_dw}。
\end{proof}

然后创建两个新的参数序列$\tilde{c}_w$和$\tilde{d}_w$，其中$\tilde{d}_w$由$d_w$序列的最大值$\bar{d}$组成，即
\begin{equation}\label{eq. PROB_tilde_d_w}
	\tilde{d}_w := [\bar{d}, \bar{d}, ... ,\bar{d}].
\end{equation}

$\tilde{c}_w$的定义分两个部分，首先使$\tilde{c}_w$的初始值与$c_w$的相同，然后使用与公式~\eqref{eq.PROB_c_w}相同的方法计算$\tilde{c}_w$其他的值，即
\begin{equation} \label{eq. PROB_tilde_c_w}
	\begin{aligned}
		\tilde{c}_0 &= c_0 = exp(-\frac{1}{2}\sum_{j=1}^{2}b_j^2)\prod_{j=1}^2(2\lambda_j)^{-\frac{1}{2}}, \\
		\tilde{c}_w &= \frac{1}{w}\sum_{r=0}^{w-1} \tilde{d}_{w-r}\tilde{c}_r, \quad w \ge 1,
	\end{aligned}
\end{equation}

$\tilde{c}_w$和$\tilde{d}_w$是原始参数序列$c_w,d_w$包络上界。由推论~\ref{cor. PROB_dw}可知，$d_w < \tilde{d}_w, \forall w \in \mathbb{R}^+$，因此$\tilde{d}_w$是$d_w$的包络上界。下面的推论~\ref{cor. PROB_recursive_c_w}给出了$\tilde{c}_w$是$c_w$的大小关系，并且推导出了$\tilde{c}_w$的一个递归计算的方法。

\begin{corollary} \label{cor. PROB_recursive_c_w}
	由式~\eqref{eq. PROB_tilde_d_w}和式~\eqref{eq. PROB_tilde_c_w}定义的序列 $\tilde{c}_w$ 的递归计算方程为：
	\begin{equation} \label{eq. PROB_c_w_recursive}
		\tilde{c}_{w+1} = \frac{\bar{d}+w}{w+1} \tilde{c}_{w}
	\end{equation}
	且$\tilde{c}_w$是式~\eqref{eq.PROB_c_w} 定义的参数序列$c_w$的包络上界，即
	$$		|c_k| \le |\tilde{c}_k|, \forall k \in \mathbb{R}^+, $$
\end{corollary}

\begin{proof}
	该推论的证明参见附录~\ref{appendix. PROB_dw}。
\end{proof}

定义$D\in \mathbb{R}^+$为恰好大于$d_w$上界$\bar{d}$的正整数，即$\bar{d} \le D < \bar{d} + 1$。根据$D$的取值，再定义一个新的序列$e_w^D$，

\begin{itemize}
	\item[1):] 当$D=1$时，	
	\begin{equation} \label{eq. PROB_e_w^1}
		e^D_w = 	\frac{c_w}{(w+1) }, \forall w \in \mathbb{R}^+,
	\end{equation}

	\item[2):] 当$D \ge 2$时，
	\begin{equation} \label{eq. new series e_k^D, D>2}
		e^D_w = \left\lbrace 
		\begin{array}{ll}
			c_w, & \forall w \le (D - 2) \\
			\frac{c_w}{\prod_{j=1}^{D} (w+2-j) }, & \forall w \ge (D - 1)
		\end{array}	\right.
	\end{equation}
\end{itemize}

则一个关于序列$e_w^D$的性质总结为如下推论：
\begin{corollary} \label{cor. PROB_e_w^D}
	假设已知正整数$D$，则序列$e_w^D$的极限是有界的，即
	$$\mathop{\lim}_{w \rightarrow \infty} |e_w^D| = \mathit{Const.} ,$$
	且$e_w^D$的最大值满足$|e_w^D| \le \tilde{c}_{D-1}, \forall w \in \mathbb{R}^+$.
\end{corollary}

\begin{proof}
	该推论的证明参见附录~\ref{appendix. PROB_e_w^D}。
\end{proof}

由以上三个推论的支撑，式~\eqref{eq. PROB_truncationError}所定义的截断误差满足如下定理，

\begin{theorem} \label{th. PROB_adaptive_main_theorem}
	给定一个二维正太向量$\vect{X} \sim \mathcal{N}(\vect{\mu},\vect{\Sigma})$，用$\vect{\lambda} = [\lambda_1,\lambda_2]^T$表示协方差矩阵$\vect{\Sigma}$的特征值，矩阵$\vect{P}$为协方差矩阵$\vect{\Sigma}$的特征向量矩阵。定义向量$\vect{b} = \vect{P}^T\vect{\Sigma}^{\frac{1}{2}}\vect{\mu} = [b_1,b_2]^T$，并且按照式~\eqref{eq.PROB_c_w}定义参数序列$c_w$，以及按照式~\eqref{eq.PROB_d_w}定义$d_w$，令$D \in \mathbb{R}^+$表示$d_w$上界的正整数。首先假设
	$$2\lambda_j >1,\forall j \in \{1,2\}, $$
	则当任意给定一个误差阈值$\delta > 0$，由式~\eqref{eq. PROB_truncationError}所定义的截断误差$\Delta p$ 可以被减少到该误差阈值之下，即
	\begin{equation*}
		\Delta p(s_{1,2}^{k+1} = 1) < \delta
	\end{equation*}
	仅当有限项近似中的累加项数目$w_F$满足如下条件时：
	\begin{equation}  \label{eq. PROB_w_F_conditions}
		\left\lbrace
		\begin{aligned}
			&w_F \ge \frac{1}{b_j^2} + \frac{1}{2\lambda_j - 1}, \forall j \in \{ 1,2 \} \\
			&w_F \ge y(\mathit{e})^2 + D - 2, \\
			& w_F \ge 3D - ln \left( \frac{\delta}{c_0y^D} \right) - 2.	 	
		\end{aligned}
		\right.		
	\end{equation}
	式中$\mathit{e}$表示自然常数。
\end{theorem}

\begin{proof}
	基于推论~\ref{cor. PROB_dw},~\ref{cor. PROB_recursive_c_w},~\ref{cor. PROB_e_w^D}，则式~\eqref{eq. PROB_truncationError}定义的截断误差可作如下展开：
	\begin{equation*}
		\begin{aligned}
			\Delta &p(s_{1,2}^{k+1}) = \sum_{w=w_F+1}^{\infty}(-1)^w c_w \frac{y^{w+1}}{(w+1)!} \\
			&\le  \sum_{w=w_F+1}^{\infty} \frac{|c_w|}{\prod_{j=1}^{D}(w+2-j)} \frac{y^{w+1}}{(w+1-D)!} \\
			&\le \sum_{w=w_F+1}^{\infty} |\tilde{e}^D_w|\frac{y^{w+1}}{(w+1-D)!} \\
			&\le \tilde{c}_{D-1} \frac{y^{w_F+2}}{(w_F+2-D)!} \left(1 + \sum_{w=1}^{\infty} \frac{y^{w}}{\prod_{j=1}^{w} (w_F + 2+ j -D) }  \right) \\
			&< \tilde{c}_{D-1} \frac{y^{w_F+2}}{(w_F+2-D)!} \left(1 + \sum_{w=1}^{\infty} \frac{y^{w}}{ (w_F -D)^w }  \right).
		\end{aligned}
	\end{equation*}
	最后一行不等式右侧的无穷累加部分可以通过巴塞尔问题来求解，即寻找合适的$w_F$使
	$$		\frac{1}{ \left( \frac{w_F-D}{y} \right)^k } \le \frac{1}{k^2},$$
	则根据巴塞尔问题可得，
	$$		\left(1 + \sum_{w=1}^{\infty} \frac{y^{w}}{ (w_F -D)^w }  \right) \le \left( 1 + \sum_{w=1}^{\infty} \frac{1}{w^2} \right)  = 1 + \frac{\pi^2}{6}. $$
	因此使得该无穷累计项满足巴塞尔问题的条件为，
	\begin{equation*}
		\begin{aligned}
			&\frac{1}{ \left( \frac{w_F-D}{y} \right)^w } \le \frac{1}{w^2}\\
			\Rightarrow\quad & \left( \frac{w_F-D}{y} \right)^w \ge w^2 \\
			\Rightarrow\quad & \ln(w_F-D) - \ln(y) \ge \frac{2\ln(w)}{w}. 	
		\end{aligned}
	\end{equation*}
	由于函数$\frac{\ln(w)}{w}$的最大值为$\frac{1}{\mathit{e}}$， 由此可得
	\begin{equation} \label{eq. PROB_w_F_Basel_condition}
		\begin{aligned}
			&\ln(w_F - D) - \ln(y) \ge \frac{2}{e} \\
			\Rightarrow \quad & w_F \ge e^{\frac{2}{e}}y + D.
		\end{aligned}
	\end{equation}
	
	\begin{table}[htbp]
		\bicaption[tab. PROB_g_w_F]{}{不同条件下$g(w_F)$低于四组域值时所需的$w_F$值}{Table$\!$}{Values of $w_F$ when $g(w_F)$ is less than a certain threshold given finite $y$ and $D$}
		\vspace{0.5em}
		\centering
		\wuhao		
		\begin{tabular}{ccccccc}
			\toprule[1.5pt]
			\multirow{2}{*}{Conditions} & \multicolumn{4}{c}{ $g(k_m) \le$ } \\ \cline{2-5}
			& $10^{-3}$ & $10^{-5}$ & $10^{-10}$ & $10^{-20}$ \\ 
			\midrule[1pt]
			y = 25($\rho$ = 5), D = 10               & 108   & 109   & 118    & 133    \\
			y = 25($\rho$ = 5), D = 20              & 139   & 141   & 149    & 163    \\
			y = 100($\rho$ = 10), D = 10              & 327  & 327  & 340   & 359   \\
			y = 100($\rho$ = 10), D = 20             & 375  & 378  & 387   & 405   \\ 
			\bottomrule[1.5pt]
		\end{tabular}
	\end{table}

	下面考虑截断误差展开的不等式的中间项，令$g(w_F) = \frac{y^{w_F+2}}{(w_F+2-D)!}w_Fw_F$，则将该序列相邻两项元素做差可得，
	$$
	\begin{aligned}
		& g(w_F+1) - g(w_F)\\ 
		&= \frac{y^{w_F+3}}{(k_m+3-D)!} - \frac{y^{w_F+2}}{(w_F+2-D)!}\\
		& = \frac{y^{w_F+2}}{(w_F+2-D)!} \left( \frac{y}{w_F+3-D} - 1\right).
	\end{aligned}
	$$
	因此，当$y=\rho^2$和$D$全都已知时，令$w_F>y+D-3$即可让$g(w_F+1)<g(w_F)$，因此$w_F$的取值越大，$g(w_F)$就越小。$g(w_F)$随$w_F$的增加而变小的趋势如表~\ref{tab. PROB_g_w_F}所示。
	
	在满足巴塞尔条件~\eqref{eq. PROB_w_F_Basel_condition}后，截断误差为
	\begin{equation} \label{eq. PROB_trunctionErrorAfterBasel}
		\begin{aligned}
			\Delta p(s_{1,2}^{k+1}) 
			&< \tilde{c}_{D-1} \left( 1+\frac{\pi^2}{6} \right)g(w_F). 
		\end{aligned}
	\end{equation}
	由推论~\ref{cor. PROB_recursive_c_w} 可知，$\tilde{c}_{D-1}$是有限值，因此经过上述的定性分析之后可得，增加$w_F$的取值可以减少截断误差。
	
	下面将定量地分析$w_F$的取值如何影响截断误差，即确定合适的$w_F$，当给定一个截断误差阈值的$\delta$时，使得$\Delta p(s_{1,2}^{k+1}) < \delta$。
	
	首先分析$\tilde{c}_{D-1}$，根据推论~\ref{cor. PROB_recursive_c_w}可知
	\begin{equation*}
		\begin{aligned}
			\tilde{c}_{w} & = \frac{w-1+\bar{d}}{w}\tilde{c}_{w-1} \\
			& = \frac{(w-1+\bar{d})(w-2+\bar{d})}{w(w-1)}\tilde{c}_{w-1} \\
			&= \frac{\prod_{i=0}^{w-1} (\bar{d}+i) }{\prod_{j=0}^{w-1} (j+1) } \tilde{c}_0. 	
		\end{aligned}
	\end{equation*}
	由于$\bar{d} < D$，$ \prod_{i=0}^{w-1} (\bar{d}+i) \le \prod_{i=0}^{w-1} (D+i) = \frac{(D+w-1)!}{(D-1)!}$，因此可知
	\begin{equation*}
		\tilde{c}_{w} \le \frac{(D+w-1)!}{w!(D-1)!} \tilde{c}_0, \forall w \in \mathbb{R}^+,
	\end{equation*}
	因此$\tilde{c}_{D-1}$的上界为
	\begin{equation*}
		\tilde{c}_{D-1} \le \frac{(2D-2)!}{(D-1)!(D-1)!} \tilde{c}_0,
	\end{equation*}	
	由引理~\ref{lem. PROB_factorial_lemma}可得
 	\begin{equation} \label{eq. PROB_maximum_tile_c_D_1}
		\begin{aligned}
			\tilde{c}_{D-1} &\le \frac{(2D-2)!}{(D-1)!(D-1)!} \tilde{c}_0 \\
			& \le \frac{e(D-1)^{2D-2}}{ \left( \frac{D-1}{e} \right)^{D-1} \left( \frac{D-1}{e} \right)^{D-1} } \tilde{c}_0 \\
			& \le e^{2D-1} \tilde{c}_0 = e^{2D-1} c_0.
		\end{aligned}
	\end{equation}	
	
	下面分析$g(w_F)$与$w_F$取值的具体关系。使用变量替换$\alpha = w_F + 2 - D$，则再次引用定理~\ref{lem. PROB_factorial_lemma}可得：
	\begin{equation} \label{eq. PROB_g(w_F)_alpha}
		\begin{aligned}
			& g(w_F) = \frac{y^{\alpha+D}}{\alpha!} \\
			& \le \frac{y^{\alpha+D}e^{\alpha}}{{\alpha}^{\alpha}} = \left( \frac{ye}{\alpha} \right)^{\alpha} y^D.
		\end{aligned}		
	\end{equation}
	将公式~\eqref{eq. PROB_g(w_F)_alpha}和公式~\eqref{eq. PROB_maximum_tile_c_D_1}代入公式~\eqref{eq. PROB_trunctionErrorAfterBasel}，则当给定截断误差的阈值$\delta$时，为了使$\Delta p(s_{1,2}^{k+1}) < \delta$，需保证
	\begin{equation*}
		\Delta p(s_{1,2}^{k+1}) \le e^{2D-1}c_0 (1+\frac{\pi^2}{6}) \left( \frac{ye}{\alpha} \right)^{\alpha} y^D < \delta.
	\end{equation*} 
	由于$1+\frac{\pi^2}{6} < \mathit{e}$，因此上式可转换为：
	\begin{equation*}
		\begin{aligned}
			e^{2D} \left( \frac{ye}{\alpha} \right)^{\alpha} < \frac{\delta f}{c_0 y^D}.
		\end{aligned}
	\end{equation*}
	
	首先令$\alpha \ge e^2y \Rightarrow w_F \ge e^2y + D - 2$,此时
	\begin{equation*}
		e^{2D} \left(\frac{ye}{\alpha} \right)^\alpha \le e^{(2D-\alpha)},
	\end{equation*}
	然后令
	\begin{equation*}
		\begin{aligned}
			& e^{(2D-\alpha)} < \frac{\delta f}{c_0y^D} \\
			\Rightarrow & \alpha > 2D - ln \left( \frac{\delta f}{c_0y^D} \right) \\
			\Rightarrow & w_F > 3D - ln \left( \frac{\delta f}{c_0y^D} \right) - 2.
		\end{aligned}
	\end{equation*}

	因此，综合巴塞尔问题对$w_F$的要求，可以总结出令$\Delta p(s_{1,2}^{k+1}) < \delta$时$w_F$的选取规则，如式~\eqref{eq. PROB_w_F_conditions}所示。至此，定理~\ref{th. PROB_adaptive_main_theorem}得证。
\end{proof}

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.6\textwidth]{figures_PROB/Trace}
	\bicaption[fig. PROB_Trace]{}{平移膨胀的协方差校正方法}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

\BiSubsection{基于平移膨胀的协方差校正}{arg2}

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwResult{ $\vect{\mu}_2,\vect{\Sigma}_2,\rho_2$ \; }
	\KwData{ $\vect{\mu}_1,\vect{\Sigma}_1,\rho_1$ \;}
	
	$\lambda_{min} \leftarrow$ 计算 $\vect{\Sigma}_1$ 的最小特征根  \;
	
	\eIf{$\lambda_{min} < \frac{1}{2}$} {
		$\beta = \frac{1}{\lambda_{min}},$
		$\vect{\Sigma}_2 = \beta \vect{\Sigma}_1,$
		$\vect{\mu}_2 = \sqrt{\beta} \vect{\mu}_1,$ 
		$\rho_2 = \sqrt{\beta}  \rho_{1}$ \;
	}
	{		$\vect{\Sigma}_2 =  \vect{\Sigma}_1,$
		$\vect{\mu}_2 =  \vect{\mu}_1$ 
		$\rho_2 =  \rho_{1}$ \;
	}
	
	
	\textbf{Return} $\vect{\mu}_2,\vect{\Sigma}_2,\rho_2$
	\AlgoBiCaption{经过修正后的有限项近似概率计算算法}{Controller for $i$-th robot at timestep $k$}
	\label{alg. Trace}
\end{algorithm}

定理~\ref{th. PROB_adaptive_main_theorem}表明，基于有限项近似的概率计算方法在选取合适的累加项个数$w_F$后，其计算结果可以任意地接近真值。然而以上的分析全部基于协方差矩阵特征值大于$\frac{1}{2}$的假设。在实际的多智能体协同任务中，智能体对未来状态的预测是实时改变的，因此协方差矩阵特征值也是实时变化的。虽然前一节的协方差近似膨胀方法在一定程度上缓解了该问题。然而该方法中多次引入了近似和假设，因此降低了整体概率预测方法的性能。本节提出一种改进的基于平移膨胀的协方差校正方法，在协方差膨胀过程中，协方差椭圆不再是原地膨胀，而是沿着圆心与坐标原点的连线作平移膨胀。因此避免了协方差近似膨胀方法中的共线假设。

该方法如图~\ref{fig. PROB_Trace}所示，图中$O_1$为膨胀之前的协方差$3\delta$置信椭圆，$O_2$为膨胀后沿着$O_1O$方向平移的结果。根据相关的几何知识易得，图中$\{O,O_1,O_2\},\{O,A,D\},\{O,C,F\},\{O,B,E\}$等四组点群，每组点群内的点都处在同一直线上。则根据几何相似可知，
\begin{equation}
	\frac{S(\widearc{ABC} \cap E_{r2})}{S(E_{r2})} =  \frac{S(\widearc{DEF} \cap E_{r1})}{S(E_{r1})}. 
\end{equation} 
因此，基于平移膨胀的协方差校正方法，在协方差矩阵膨胀后不需要做任何近似和假设，即可以使用膨胀后新的且满足定理~\ref{th. PROB_adaptive_main_theorem}的增广协方差矩阵进行有限项近似的概率计算，此时的等效半径$\rho_{2}$将大于最初的工作半径$\rho_{1}$。

本节所提出的协方差平移膨胀和前一节所提出的近似屏障方法都是基于高斯噪声假设。只有在高斯噪声条件下才可以使用椭圆置信区间来表示协方差矩阵。因为在高斯模型下理论上将会有超过$99.7\%$的抽样点落在置信椭圆内。

基于平移膨胀的方法也可以被反向使用，即当通信和观测半径$\rho_{1}$比较大，但是协方差矩阵的分布较为分散时，可以反向地使用该膨胀方法，把一个分散的协方差矩阵缩小但是依然符合定理~\ref{th. PROB_adaptive_main_theorem}对协方差矩阵的要求，但是可以换取一个较小的等效半径$\rho_{2}$。该方法的好处是由于$y=\rho^2$在有限项近似的计算公式中被频繁使用，一个较大的工作半径$\rho_{1}$会积累数值计算误差，从而导致最终的结果偏离正确概率。因此较小的$\rho_{2}$将减轻此类问题。

下面采用与协方差近似膨胀相同的方法计算膨胀系数$\beta$，即确保协方差矩阵$\vect{\Sigma}$的最小特征值不小于1。下面将该算法的伪码总结为算法~\ref{alg. Trace}。


\begin{algorithm*}[htbp]
	\SetAlgoLined
	\KwResult{ 	$\qquad p(s_{1,2}^{k+1}=1)$  \;}
	\KwData{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}, \rho_{1}, \delta$ \;}
	
	$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
	$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
	
	$[\rho_{3\delta}^{min} \rho_{3\delta}^{max}] \leftarrow 3\delta$ 置信域计算( $\vect{\Sigma}_{\Delta}$ )
	
	\If{$\rho < \rho_{3\delta}^{min}$}{
		\textbf{Return} $p(Y<\rho^2) = 0$ \;
	} 
	\If{$\rho > \rho_{3\delta}^{max}$}{
		\textbf{Return} $p(Y<\rho^2) = 1$ \;
	} 
	$\left[ \vect{\mu},\vect{\Sigma},\rho_2 \right] \leftarrow \text{调用算法~\ref{alg. Trace}对协方差矩阵进行校正}(\vect{\mu}_{\Delta},\vect{\Sigma}_{\Delta},\rho_1)$ \;
	
	$\vect{\lambda} = [ \lambda_1, \lambda_2 ]^T \leftarrow$ 计算协方差矩阵 $\vect{\Sigma}$的特征值 \;
	
	$\vect{P} \leftarrow$ 计算协方差矩阵$\vect{\Sigma}$的特征向量矩阵 \;
	$\vect{b} = [b_1,b_2]^T = P^T\vect{\Sigma}^{-\frac{1}{2}}\vect{\mu}$\;
	
	$y = \rho^2 $, $w_F = 50$ \tcc*{赋予$w_F$初始值} 
	
	\For{$r=0;r\le w_F;r=r+1$}{
		$d_r \leftarrow $ 调用公式~\eqref{eq.PROB_d_w}, $\forall r = 1,...,w_F$\;
		$c_r \leftarrow $ 使用 $d_r$ 并调用公式~\eqref{eq.PROB_c_w} $\forall r = 0,...,w_F$\;
		$y_r = (-1)^r \dfrac{y^{r+1}}{\Gamma(r+2)}, r = 0,...,w_F $  \;
		\If{$d_r$达到最大值}{
			$D = d_r$ \tcc*{该条件语句仅进入并执行一次}
			$w_{F1} = \frac{1}{b_1^2} + \frac{1}{2\lambda_1 - 1}$,$w_{F2} = \frac{1}{b_2^2} + \frac{1}{2\lambda_2 - 1}$ \;
			$w_{F3} = e^2y + D-1$, $w_{F4} = 3D - ln \left( \frac{\delta}{c_0y^D} \right) - 2$ \;
			
			$w_F = \text{max}(k_{m1},k_{m2},k_{m3},k_{m4},k_m)$ \;
		}
	}
	$ p(s_{1,2}^{k+1}=1) \leftarrow $ 基于 $c_r, y_r$调用公式~\eqref{eq. PROB_CDP_finite}预测连通概率 \;
	
	\textbf{Return} $p(s_{1,2}^{k+1}=1)$
	\AlgoBiCaption{基于自适应参数选取的有限项近似概率改进算法}{Controller for $i$-th robot at timestep $k$}
	\label{alg. PROB_APSE}
\end{algorithm*}

\BiSubsection{基于自适应参数选择的连通概率预测算法}{arg2}

前两节给出了两种有限项近似的改进策略。首先是通过严格的理论推导，得到了有限项近似方法累加项数目$w_F$的自适应选取策略。然后是针对自适应策略中存在的对协方差矩阵的约束，设计了一种改进的协方差膨胀方法。相比于前一节提出的协方差近似膨胀方法，改进的膨胀策略避免了近似共线假设，进而提升了有限项近似在协方差矩阵太小或者近似奇异情况下的计算性能。基于以上分析，将基于自适应参数选取的改进概率预测算法的伪码总结为算法~\ref{alg. PROB_APSE}。

按照从上到下的顺序，算法~\ref{alg. PROB_APSE}首先需要两个节点对未来$k+1$时刻的预测状态$\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$以及节点上装配的相对通信和测量设备的工作距离$\rho$和期望的截断误差精度$\delta$作为算法输入。最终返回两个节点在未来时刻存在通信和相对测量连通的概率，即$p(s_{1,2}^{k+1})$。算法的$3-9$行是前一种有限项近似的概率计算方法中的$3\delta$区间校正。保留该修正过程的目的是为了在节点之间的相对距离远大于或远小于设备的工作半径时，直接判定连通度的存在与否以减少算法的计算复杂度。然后第$10$行调用了改进的平移膨胀协方差校正算法，以确保节点之间相对距离分布的协方差矩阵符合定理~\ref{th. PROB_adaptive_main_theorem}的要求。$11-13$行是有限项近似公式~\eqref{eq. PROB_CDP_finite}的前期变量处理过程。$14-26$行执行有限项近似的概率计算过程，其中$19-24$行为有限项近似中累加项数目$w_F$的自适应选择方法，对应于定理~\ref{th. PROB_adaptive_main_theorem}中的选取条件。在自适应地选取$w_F$的值之前，算法为$w_F$选定了一个初始值，即$14$行中的$w_F=50$。该初始化过程实际上不是必需的，仅是为了初始情况下方便程序对内存进行分配。另外，该算法假设在给定的初始累加项数目里（此算法是50项），$d_w$能够寻找到该序列的最大值，因此可以进入$19-24$行的条件语句。由于最终的$w_F$的取值是将定理~\ref{th. PROB_adaptive_main_theorem}中所有的候选值与初始值作比较并取最大值，因此有限项近似的计算公式中将至少累加50项。

\BiSubsection{数值仿真验证}{arg2}

本节测试算法~\ref{alg. PROB_APSE}对节点之间连通存在概率预测的性能表现，此处考虑两类性能指标，即算法对连通概率预测的准确性和算法的数值运算复杂度。为了验证算法的有效性，此处的仿真验证将算法~\ref{alg. PROB_APSE}与采取不同抽样次数的随机抽样方法作对比。随机抽样方法是一种通过拟合样本分布，经过重复抽样，统计事件发生的频率来模拟概率的方法。将使用随机抽样方法预测连通概率的过程总结为算法~\ref{alg. RandomSampleDistributionAlg}。

\begin{algorithm}[htbp]
	\SetAlgoLined
	\KwResult{ $p(s_{1,2}^{k+1}=1)$  }
	\KwData{ $\vect{\tilde{p}}_1^{k+1}, \vect{\tilde{p}}_2^{k+1},\vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}}, \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}},\rho dim$ \;}

	$\Delta \vect{p}_{12}^{k+1} \sim \mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$\;
	$\qquad \vect{\mu}_{\Delta} = \vect{\tilde{p}}_1^{k+1} - \vect{\tilde{p}}_2^{k+1}$; $\vect{\Sigma}_{\Delta} = \vect{\Sigma}_{\vect{\tilde{p}}_1^{k+1}} + \vect{\Sigma}_{\vect{\tilde{p}}_2^{k+1}}$ \;
	
	Count = 0 \;
	\While{$k \le dim$}{
		$\vect{\mu}$ $\leftarrow$ 根据$\mathcal{N}(\vect{\mu}_{\Delta}, \vect{\Sigma}_{\Delta})$的分布进行随机抽样 \;
		\If{ $||\vect{\mu}||_2$ $ \le \rho$}
		{Count = Count + 1 \;}
	}
	$p(s_{1,2}^{k+1}=1) = \frac{Count}{dim}$\;
	
	\textbf{Return} $p(s_{1,2}^{k+1}=1)$
	\AlgoBiCaption{随机抽样算法}{Controller for $i$-th robot at timestep $k$}
	\label{alg. RandomSampleDistributionAlg}
\end{algorithm}

算法~\ref{alg. PROB_APSE}和算法~\ref{alg. RandomSampleDistributionAlg}的仿真实现是基于 Windows 系统下的 MATLAB R2019a 软件，运行于i7-9750H@2.6GHz 处理器和8GB内存的计算平台。由于算法~\ref{alg. PROB_APSE}在具体的程序实现上是运动MATLAB的\textit{for-loop}语句，而算法~\ref{alg. RandomSampleDistributionAlg}的实现是采用MATLAB软件的内置函数\textit{mvrnd}函数，因此算法~\ref{alg. PROB_APSE}在数值计算效率上存在算法层面和优化层面的提升空间。

下面将仿真的具体设置介绍如下。测试两个节点$A$和$B$，它们的真实位置固定为$\vect{p_A} = [0.5,1]^T$ 以及 $\vect{p}_B=[2,2.5]^T$，因此节点$A，B$之间的真实距离为$||\vect{p}_A - \vect{p}_B||_2 \approx 2.21$。然后随机生成两个正定对称矩阵$\vect{P}_A>0,\vect{P}_B>0$作为对$A,B$两个节点估计位置分布的协方差矩阵。对任意的一组协方差分布，设备的工作半径共测试$\rho$为$\rho \in \vect{S}_{\rho}:=\{ 0.1:0.1:6 \}$集合中的60个数值，以此验证算法对连通存在的确定存在状态、确定不存在状态以及临界状态的表现性能。

\begin{figure}[htbp]
	\centering
	\includegraphics[width = 0.9\textwidth]{figures_PROB/32RealDistributionAndErrorsSingleCase}
	\bicaption[fig. PROB_RealDistributionAndErrorsSingleCase]{}{典型案列的真实分布}{Fig.$\!$}{Scenairo of active and cooperative localization}\vspace{0em}
\end{figure}

\begin{figure}[htbp]
	\centering	
	\subfigure{\label{subfig. discretizationArea}}\addtocounter{subfigure}{-2}
	\subfigure[leftfig]{\subfigure[RMSE分布]{\includegraphics[width = 0.45\textwidth]{figures_PROB/32RMSESingleCase}}}	
	\null\hfill
	\subfigure{\label{subfig. discretizationBoundary}}\addtocounter{subfigure}{-2}
	\subfigure[leftfig]{\subfigure[数值计算复杂度]{\includegraphics[width = 0.45\textwidth]{figures_PROB/32RunningTimeSingleCase}}}
	
	\bicaption[fig. PROB_performanceRepresentativeCase]{}{典型案例下算法的性能表现与对比}{Fig.$\!$}{Gradient of the objective on the coverage boundary}
	\vspace{-0em}
\end{figure}

然后对于任意的$\rho \in \vect{S}_{\rho}$，将$\vect{p}_A,\vect{p}_B,\vect{P}_A,\vect{P}_B$和$\rho$作为算法~\ref{alg. PROB_APSE}以及算法~\ref{alg. RandomSampleDistributionAlg}的输入。对于算法~\ref{alg. PROB_APSE}最后的待输入参数，截断误差的期望阈值$\delta$，此处设置为固定值$10^{-10}$。随机抽样算法~\ref{alg. RandomSampleDistributionAlg}共测试六组不同的抽样值作为对照组，即$dim=[10,10^2,10^3,10^4,10^5,10^6]$。概率的真实值使用$dim=10^7$的算法~\ref{alg. RandomSampleDistributionAlg}得出的频率值来近似。


首先测试一组典型案列，随机生成两个协方差矩阵为
$$\vect{P}_A = \left[ \begin{matrix}
	0.371629432909261&-0.100755606016948 \\
	-0.100755606016948&	0.263717689138266
\end{matrix} \right],$$
以及
$$\vect{P}_B = \left[ \begin{matrix}
	0.196735055559263&0.370216978853931 \\
	0.370216978853931&0.720201443148860
\end{matrix} \right].$$
因此节点$A$和$B$之间的相对距离分布为$\Delta \vect{p} \sim \mathcal{N}(\vect{p_A} - \vect{p}_B, \vect{P}_A+\vect{P}_B)$。通过计算可知，协方差$\vect{P}_A+\vect{P}_B$的最小特征值为 $0.435875962912945<0.5$，因此，该典型案例在调用算法~\ref{alg. PROB_APSE}时将激活算法~\ref{alg. Trace}进行协方差膨胀的校正过程。

典型案例中两个节点的分布如图~\ref{fig. PROB_RealDistributionAndErrorsSingleCase} A 所示，连个节点的相对距离分布如图中 B 所示。子图 C 展示了实验组和6个对照组概率预测结果相对于真实概率的误差随设备工作半径$\rho$变化的分布。对实验组和6个对照组对$\vect{S}_{\rho}$中60个半径的仿真结果作处理，可分别获得相对于真实概率的概率预测误差的RMSE值以及每组方法的平均运算时间（包括真实概率的仿真组），此两种性能表现的对比参见图~\ref{fig. PROB_performanceRepresentativeCase}。

\begin{figure}[htbp]
	\centering	
	\subfigure{\label{subfig. discretizationArea}}\addtocounter{subfigure}{-2}
	\subfigure[leftfig]{\subfigure[概率误差RMSE分布]{\includegraphics[width = 0.45\textwidth]{figures_PROB/81RMSEALLCases}}}	
	\null\hfill
	\subfigure{\label{subfig. discretizationBoundary}}\addtocounter{subfigure}{-2}
	\subfigure[leftfig]{\subfigure[数值计算复杂度分布]{\includegraphics[width = 0.45\textwidth]{figures_PROB/81RunningTimeALLCases}}}
	
	\bicaption[fig. PROB_performance200Distribution]{}{200次仿真中算法的概率预测表现性能分布}{Fig.$\!$}{Gradient of the objective on the coverage boundary}
	\vspace{-0em}
\end{figure}

然后将上述典型案例的仿真过程重复200次，每次随机生成节点的位置状态协方差矩阵$\vect{P}_A,\vect{P}_B$。计算每一次仿真对60个设备工作半径$\rho$的平均概率预测误差的RMSE值以及包括真实概率组在内的每组方法的平均运算时间。此两种性能表现在200次的重复实验中的分布如图~\ref{fig. PROB_performance200Distribution}。

由仿真结果可得，本节设计的基于自适应参数选取的有限项近似的概率预测算法~\ref{alg. PROB_APSE}在概率预测精度上优于所有的对照组。其数值计算复杂度上介于使用$10^{3}$-$10^4$个抽样的随机抽样算法的计算复杂度之间。

\BiSection{小结}{Conclusion}

本章研究了多智能体系统如何预测在未来时刻系统内部的通信观测连通概率问题。该问题源自系统对未来时刻系统行为进行预测时存在，系统的未来时刻状态在当前时刻是不确定、随机的。本章基于系统对未来时刻状态值的预测分布，针对通信和观测过程的圆盘模型，设计了一类可以预测任意两个节点在未来时刻存在通信和观测链概率的方法。该方法基于正太二次型随机变量的多项式展开定理，并通过使用有限项近似来实现数值计算的可行性。通过严格的理论推导和实验仿真验证，设计了一系列的修正方法和自适应的参数选择策略，用以克服因有限项近似引入的截断误差的影响。本章最终设计的自适应参数选取的概率预测算法具有如下特性：
\begin{itemize}
	\item[1):] 严格的理论推导和数值仿真证明了有限项近似所引入的截断误差可以通过选取合适的多项式累加数而任意减少。
	\item[2):] 该算法最终的计算复杂度介于使用$10^3$-$10^4$次抽样的随机抽样算法的复杂度之间。
\end{itemize}

本章设计的连通概率预测方法为后文对多智能体系统主动协同定位方法的研究提供了支撑。

然而本章研究的方法依然存在不足。最重要的缺点是数值计算平台的计算精度问题。因为在有限项近似的计算公式中，变量$y=\rho^2$的高次幂运算被频繁地使用。然而由于最终的概率结果仅限制在概率的值域 $[0,1]$ 之间，因此最终的计算结果很容易受到计算平台数值精度的影响。未来的研究方向可标准卡方分布来对相对距离的二次型分布进行建模，从而可以使用 Patnaik \citeup{patnaik1949non} 近似或者 Pearson \citeup{imhof1961computing} 的近似方法，通过查表来获取最终的概率值。该方法无需大量的幂次运算，但是最终的精度或将受到查表精度的影响。

\BiSection{附录}{}

\BiSubsection{证明推论~\ref{cor. PROB_dw}}{arg2} \label{appendix. PROB_dw}
\begin{proof}
	根据式~\eqref{eq.PROB_d_w}中对$d_w$的定义，将该序列的两个相邻项作减法可得：
	\begin{equation*}
		\begin{aligned}
			&d_{w+1} - d_w  \\
			&= \frac{(2\lambda_j)^{-w}}{2}  \sum_{j=1}^2 \left\lbrace \left[ 1-(w+1)b_j^2 \right] (2\lambda_j)^{-1} - (1-wb_j^2) \right\rbrace  \\
			&= \frac{(2\lambda_j)^{-w}}{2}  \sum_{j=1}^2 \left\lbrace  wb_j^2\left[ 1 - (2\lambda_j)^{-1} \right] + (1-b_j^2)(2\lambda_j)^{-1} - 1\right\rbrace.
		\end{aligned}
	\end{equation*}
	由于假设$\left[ 1 - (2\lambda_j)^{-1} \right]>0$，因此如果
	\begin{equation*} 
		\begin{aligned}
			w &\ge \frac{1 + (b_j^2 - 1)(2\lambda_j)^{-1}}{ b_j^2\left[ 1 - (2\lambda_j)^{-1} \right]}, \\
			& \ge \frac{1}{b_j^2} + \frac{1}{2\lambda_j - 1}
		\end{aligned},  \forall j \in \{ 1,2 \},
	\end{equation*}
	$d_w$序列相邻两项的差值有
	$$d_{w+1} - d_{w} > 0.$$
	即$d_w$序列在$w$大于某个阈值后将持续增加。
	
	另外，从$d_w$的定义可知，当
	\begin{equation}
		w \ge \frac{1}{b_j^2},  \forall j \in \{ 1,2 \}. 
	\end{equation}
	时，有$d_w \le 0, \forall k \in \mathbb{R}^+$。因此可知，序列$d_w$在$w$足够大时，随着$w$的变大而增加的负数。因此该序列的极限必然是无限趋近于0。
	
	又因为$d_w$对任意的$w$都有定义，$d_1\neq \infty$，因此必然存在某一元素是整个序列的最大值$\bar{d}$，使得$|d_w|\le \bar{d}, \forall w \in \mathbb{R}^+$。
	
	至此，推论~\ref{cor. PROB_dw}得证。
\end{proof}


\BiSubsection{证明推论~\ref{cor. PROB_recursive_c_w}}{arg2} \label{appendix. PROB_recursive_c_w}
\begin{proof}
	首先将$\tilde{c}_w$序列的两个相邻项$w$和$w+1$进行展开，
	 \begin{equation} \label{eq. PROBappen_expan_c_w}
		w\tilde{c}_w = \left( \tilde{c}_0\tilde{d}_w + \tilde{c}_1\tilde{d}_{w-1} + ... + \tilde{c}_{w-1}\tilde{d}_1\right),
	\end{equation}
	以及
	\begin{equation} \label{eq. PROBappen_expan_c_w+1}
		(w+1) \tilde{c}_{w+1} = \left( \tilde{c}_0\tilde{d}_{w+1} + \tilde{c}_1\tilde{d}_{w} + ... + \tilde{c}_{w}\tilde{d}_1\right).
	\end{equation}
	则将式~\eqref{eq. PROBappen_expan_c_w}的左右两边同时添加一个$\tilde{c}_w$可得，
	 \begin{equation} \label{eq. PROBappen_expan_c_w+1_add_onemoreterm}
		(w+1)\tilde{c}_w = \left( \tilde{c}_0\tilde{d}_w + \tilde{c}_1\tilde{d}_{w-1} + ... + \tilde{c}_{w-1}\tilde{d}_1 + \tilde{c}_w \right).
	\end{equation}
	用式~\eqref{eq. PROBappen_expan_c_w+1} 减去式~\eqref{eq. PROBappen_expan_c_w+1_add_onemoreterm}可得，
	 \begin{equation*}
		\begin{aligned}
			&(w+1)\left( \tilde{c}_{w+1} - \tilde{c}_{w} \right) \\
			&= \tilde{c}_0 (\tilde{d}_{w+1} - \tilde{d}_{w} ) + ... + \tilde{c}_{w-1} (\tilde{d}_{2} - \tilde{d}_{1} )  +\tilde{c}_w (\tilde{d}_{1} - 1 ).  
		\end{aligned}		
	\end{equation*}
	由于$\tilde{d}_{w+1} = \tilde{d}_{w} = ... = \tilde{d}_{1} = \bar{d}$，因此上式转变为：
	 \begin{equation*} 
		\begin{aligned}
			(w+1)&\left( \tilde{c}_{w+1} - \tilde{c}_{w} \right) = (\bar{d} - 1)\tilde{c}_w \\
			\Rightarrow \quad & \tilde{c}_{w+1} = \frac{\bar{d}+w}{w+1} \tilde{c}_{w}. 
		\end{aligned}		
	\end{equation*}
	另外，因为$\tilde{c}_0 = c_0 > 0$，因此$\tilde{c}_w > 0, \forall k \in \mathbb{R}^+$。
	
	下面使用数学归纳法证明$\tilde{c}_w$是序列$c_w$的包络上界。
	
	首先当$i=0$时， $c_0 = \tilde{c}_0 \Rightarrow |c_0| \le |\tilde{c}_0|$。 
	 
	 然后当$i=1$时， $|c_1| = |c_0d_1| \le |c_0| |d^u| = |\tilde{c}_1|.$
	 
	 假设当$i=w-1$时，满足$|c_{w-1}| \le |\tilde{c}_{w-1}|$。 则对于$i=w$，有
	  \begin{equation*} 
	 	\begin{aligned}
	 		|c_w| &= |\frac{1}{w} \left( c_0d_w + c_1d_{w-1} + ... + c_{w-1}d_1 \right)| \\
	 		& \le \frac{1}{w} \left( |c_0||d_w| + |c_1||d_{w-1}| + ... + |c_{w-1}||d_1| \right)\\
	 		& \le \frac{1}{w} \left( |\tilde{c}_0||\bar{d}| + |\tilde{c}_1||\bar{d}| + ... + |\tilde{c}_{w-1}||\bar{d}| \right) \\
	 		& = |\tilde{c}_w|.
	 	\end{aligned}	
	 \end{equation*}
 	因此对任意的$w$,都满足$|c_w| \le |\tilde{c}_w|$。
 	
 	至此，推论~\ref{cor. PROB_recursive_c_w}得证。
\end{proof}



\BiSubsection{证明推论~\ref{cor. PROB_e_w^D}}{arg2} \label{appendix. PROB_e_w^D}

\begin{proof}
	证明新序列$e_w^D$在该推论中各种性质的思路是，使用$c_w$的上界序列$\tilde{c}_w$，依照$e_w^D$序列的定义方式，用同样的定义方法使用$\tilde{c}_w$合成一个$\tilde{e}_w^D$。由于$\tilde{c}_w$是$c_w$的上界序列，该性质可以拓展到具有相同定义方式的序列$\tilde{e}_w^D$和$e_w^D$上，即
	$$|e_k^D| \le |\tilde{e}_k^D|,\forall k \in \mathbb{R}^+.$$	
	
	则根据推论~\ref{cor. PROB_recursive_c_w}中$\tilde{c}_w$的递归计算式，当$w>D-1$时，有如下关系式的成立，
	\begin{equation*} \label{eq. tilde_e_k^D and tilde_e_{k+1}^{D}}
		\tilde{e}_{w+1}^D = \frac{\bar{d}+w}{w+1} \frac{\prod_{j=1}^{D}(w+2-j)}{\prod_{j=1}^{D}(w+3-j)}\tilde{e}_{k}^D.		
	\end{equation*}

	因此当$D=1$，上式可简化为
	\begin{equation*}
		\begin{aligned}
			\tilde{e}_{w+1}^D &= \frac{\bar{d}+w}{w+2} \tilde{e}_{w}^D \\
			& < \tilde{e}_{w}^1, \quad \forall w \ge 0.
		\end{aligned}		
	\end{equation*}
	因此当$D=1$时，序列$\tilde{e}_w^D$在随着$w$的增加而逐渐减少。此时$\tilde{e}_w^D$的最大值为$\tilde{e}_0^D = \tilde{c}_0$。
	
	当$D \ge 2$时，根据式~\eqref{eq. new series e_k^D, D>2}的定义方式，可做如下展开，	
	\begin{equation}
		\begin{aligned}
			\tilde{e}_{w+1}^D &= \frac{\bar{d}+w}{w+1} \frac{(w+1)(w+2-D)\prod_{j=2}^{D-1}(w+2-j)}{(w+2)(w+1)\prod_{j=3}^{D}(w+3-j)}\tilde{e}_{w}^D  \\
			& \le \frac{ (\bar{d}+w)(w+2-D) }{ (w+2)(w+1) } \frac{\prod_{j=2}^{D-1}(w+2-j)}{\prod_{j=2}^{D-1}(w+2-j)}\tilde{e}_{w}^D \\
			& = \frac{ (\bar{d}+w)(w+2-D) }{ (w+2)(w+1) }\tilde{e}_{w}^D\quad \forall w \ge D-1.
		\end{aligned}		
	\end{equation}
	将系数的分子分母做差可得
	$$
	\begin{aligned}
		&(\bar{d}+w)(w+2-D) - (w+2)(w+1) \\
		&= (\bar{d}-D-1)w +\bar{d}(2-D)-2.
	\end{aligned}
	$$
	由于$0 \le \bar{d} \le D$并且$D \ge 2$，因此$\bar{d} -D -1 < 0$，以及$\bar{d}(2-D) \le 0$，因此可知
	\begin{equation} \label{eq. tilde_e_k^D if D>=2}
		\tilde{e}_{w+1}^D < \tilde{e}_{w}^D, \quad \forall w \ge D-1.
	\end{equation}
	由此可以总结出$\tilde{e}_w^D$的是有界的且它的极限是有限值，即
	$$\mathop{\lim}_{k \rightarrow \infty} \tilde{e}_{k}^D = \mathit{Const.}$$
	此时当$w\le D-2$，根据式~\eqref{eq. new series e_k^D, D>2}中的定义，$\tilde{e}_w^D$的值是随着$w$的增加而递增的，因此在$0 \ge w \le D-2$的区间内，$\tilde{e}_w^D$在$w=D-2$取得最大值。又当$w \ge D-1$时，由上述分析可知，$\tilde{e}_w^D$在$w=D-1$处取得最大值。又$\tilde{e}_{D-2}^D=\tilde{c}_{D-2} < \tilde{c}_{D-1}$，并且$\tilde{e}_{D-1}^D=\frac{\tilde{c}_{D-1}}{D!} \le \tilde{c}_{D-1}$，因此可以总结出，$\tilde{e}_w^D < \tilde{c}_{D-1}, \forall w \in \mathbb{R}^+$。
	
	
	又由于$\tilde{e}_w^D$ 是 $e_w^D$的上界包络，因此$|e_w^D|$的极限也是有限值，且$|e_w^D|<\tilde{c}_{D-1}$。
	
	至此，推论~\ref{cor. PROB_e_w^D}得证。

\end{proof}